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Abstract 

We present a new method, in the following called MMM2D, to accurately calculate 
the electrostatic energy and forces on charges being distributed in a two dimensional 
periodic array of finite thickness. It is not based on an Ewald summation method and 
as such does not require any fine-tuning of an Ewald parameter for convergence. We 
transform the Coulomb sum via a convergence factor into a series of fast decaying 
functions which can be easily evaluated. Rigorous error bounds for the energies 
and the forces are derived and numerically verified. Already for small systems our 
method is much faster than the traditional 2D-Ewald methods, but for large systems 
it is clearly superior because its time demand scales like 0(iV 5 / 3 ) with the number 
N of charges considered. Moreover it shows a rapid convergence, is very precise and 
easy to handle. 
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1 Introduction 



The calculation of Coulomb interactions are time consuming due to their long 
range nature. In principle, each charge qi at position pi interacts with all oth- 
ers, leading to a computational effort of 0(N 2 ) already within the central 
simulation box. For many physical investigations one wants to simulate bulk 
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properties and therefore introduces periodic boundary conditions to avoid sur- 
face effects. The Coulomb energy E then has to be computed as a sum over 
all periodic images, 



\Pi - Pj + n\ 



(1) 



raeZ 3 i,j 



where for simplicity we assume a cubic simulation box of length A. The prime 
denotes that for n = the term % — j has to be omitted. This sum is merely 
conditionally convergent, meaning that its value depends on the summation 
order. Normally a spherical limit is used, i. e. the vectors n are added in 
the order of increasing length. Then E can be computed for example via the 
traditional Ewald summation method [1]. The basic idea of this method is 
to split the original sum via a simple transformation involving the splitting 
parameter a into two exponentially convergent parts, where the first one is 
short ranged and evaluated in real space, the other one is long ranged and can 
be analytically Fourier transformed and evaluated in Fourier space. For any 
choice of the Ewald parameter a and no truncation in the sums the formula 
yields the exact result. In practice one cuts off the infinite sums at some finite 
values and obtains E to a user controlled accuracy, which is possible by using 
error estimates for the cut-offs [2]. 

Another way of computing Eq. (1) is via a convergence factor 



This approach is used in a method due to Sperb [3] called MMM. 

In Refs.[4,5] it was shown that Eqs. (1) and (2) differ by a term proportional to 
the dipole moment for 3D systems. Sperb et al. have developed a factorization 
approach which yields an 0(N log N) algorithm [6], comparable in speed to 
the FFT mesh based Ewald algorithms [7, 8]. An different convergence factor 
was used by Lekner [9] to efficiently sum up the 3D Coulomb sum. But this 
approach still leads to an 0(N 2 ) algorithm. 

For thin polyelectrolytes films, membrane interactions, or generally confined 
charged systems, one is interested in summations where only 2 dimensions are 
periodically replicated and the third one (here always the z-direction) is of 
finite thickness h (2D + h geometry). For this geometry Ewald based formulas 
are only slowly convergent, have mostly 0(N 2 ) scalings and no "a priori" error 
estimates exist [10]. In Ref.[10] a comparison was made with respect to speed 
and accuracy of the three most popular versions, the traditional 2<i-Ewald 
method developed by Heyes, Barber, and Clarke (HBC) [11], the one due to 



= lim \ ^ 

/3^0 2 ^ 
nSZ 3 13 



qiqj exp (-f3\pj - pj + wAj) 
\pi - pj +n\\ 



(2) 
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Hautmann-Klein (HK) [12], and the approach attributed to Nijboer and de 
Wette (NdW) [13,14]. Another variant is due to Smith [15], which combines 
HBC and NdW, resulting in a good scaling, but still slow convergence. Re- 
cently two more Ewald-type methods were proposed with different ways of 
approximating the Fourier space sum of the HBC-formula[16,17]. 

The aim of this paper is to develop a new method, which we call MMM2D. 
A brief description of the key features of this algorithm was already given 
in Ref.[18]. MMM2D will be using the same ideas leading to MMM and 
is therefore derived from an equation like Eq. (2) using a convergence fac- 
tor. We will show that for two dimensional systems this approach yields the 
same results as the spherical limit, i. e. E — E. In principle, already the 3D 
formula developed by Lekner [19] could be used for 2D + h Systems. It em- 
ploys modified Bessel functions, but is only useful for particles separated in 
the ^-direction, and is only calculable pairwise. Sperb [20] derived a different 
formula, which is also applicable for particles lying in the same z-plane. A 
combination of theses two formula has been used recently in MD simulations, 
because it is easy to implement [21] for small number of particles, but suf- 
fers the drawback of being an 0(N 2 ) method, and therefore the system sizes 
which can be investigated are limited. Moreover these methods use different 
convergence factors, and up to now no investigations have been performed to 
see whether both methods lead to the same result as the one obtained by using 
a spherical limit. 

To precisely define the problem we consider systems of N particles with charges 
qi G K. and pairwise different coordinates pi = (xi,yi,Zi) T e B , i = 1, . . . , N, 
where 



B n 




2 ' 2 



xRcl 3 



is the simulation box. Furthermore it is assumed that the system is charge 
neutral, i. e. X^Li Qi = 0. 

Now let n k i := (k\ x , l\ y , 0) T for fc, I 6 Z. The expressions to be calculated are 
the interaction energies of particle i with all other charges, defined as 

oo N 

^ = E £ E' i m \ i i = i...N. (3) 

k 2 +l 2 =S 

As usual the prime ' on the inner sum indicates that the term j = i for 
k — I — is omitted. 

In Sec. 2 we will develop the formulas used by MMM2D for calculating the 
energy and the forces. In Sec. 3 we develop error formulas for the energy 
and force calculations, which can be easily calculated prior to a simulation. 



3 



Then, in Sec. 4 we give the details of our implementation, which is based on a 
factorization approach resulting in an 0(N 5 ^ 3 ) scaling. In Sec. 5 we apply our 
algorithm to several test cases and prove its applicability and efficiency. We will 
also demonstrate that it is much faster than the 2oKEwald method, and is thus 
a superior alternative to all aforementioned methods. We will end with some 
conclusions, and show for the interested reader in the Appendix A that for 
two dimensional systems E = E, i. e. that the summation using a convergence 
factor as in Eq. (2) yields exactly the same result as the summation using a 
spherical limit as the 2<i-Ewald method, there is no dipolar correction term 
needed. 



2 The MMM2D Method 



The transformations given here are similar to the transformations used in 
R. Strebel's MMM [22,6] for three-dimensional periodic systems. The ideas 
used in the following are similar to the ideas used to develop the formulas for 
MMM. Basically, two different formulas are used to calculate the energy and 
forces. For the interaction between particles that are well separated, a formula 
that can be calculated in linear time is used. For the interaction between neigh- 
boring particles a second formula is used. This leads to an 0(N 7 ^) algorithm. 
Up to here the same ideas are used for MMM2D. But for the three dimen- 
sional case some further tricks are used to achieve the scaling of 0(N log N). 
These tricks involve the symmetry of the coordinates which is not true for two 
dimensional systems. Therefore this better scaling cannot be carried over to 
the two dimensional case. 

The formulas needed for MMM2D are derived in the next subsections in the 
same way as has been done for MMM [22]. We first rewrite the sum (3) using 
a convergence factor: 

N , q.q. e -P\Pij+n k i\ n 

E * = i im n E E M —z i = ^^QiQjM^yij^i) ( 4 ) 

^° *,jez j=i mi + n ki I ^° j= i 
where p i3 - = {xij, yij, %) =Pi- Pj, 

{x,y,z) = (0,0,0) 



and 



e -@ r ki 

<l>p{x,y,z)= J2 — — , (6) 
(fc.O^Co.o) '« 
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where we abbreviated 



r k i = 



^(x + kX x ) 2 + (y + lX y ) 2 + z\ 



r k = r k0 , r = r . 



The fact that the rhs. of Eq. (4) is equal to Eq. (3) is non-trivial; see ap- 
pendix A for a proof. 



2. 1 Transformation of <pp for z ^ 

First we concentrate on developing an absolutely and rapidly converging for- 
mula for (j>p. Then we can easily form the limit (3 — > and obtain a formula 
for 0. For z ^ and (5 > the sum in <pp is an absolutely convergent sum 
of Schwartz class functions. Therefore for 5 > and x G R, we can apply the 
Poisson formula: 



where T denotes the Fourier transformation. Furthermore we will be using 
the formulas 



which are valid for a,z6l and can be found, for example, in [23]. Ko is called 
the modified Bessel function of order 0. For properties of the Bessel functions, 
see [24]. 






Setting 



P Pq = yJ{P + (2nu x pf + (2vr M ,g)2 , 



^J(3 2 + (2ttu xP ) 2 and (3 q = J(P + {2^u y q) 2 , 




l\y) 2 + Z 2 , p = po, U x = — and u y = 



1 



'X 



Xy 
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we obtain after two Fourier transformations 



e -l3r k{ I e -/3r kt 



Mw)= E — -= EE , 

k,lei kl lei \kez ' kl 



p fipq 1 2 1 

= 2nu x u y — R e 2niu * px e 2niu » qy 

p,qez Ppi 

Expanding the term for p = q = we find e ^ N = — \z\ + O n (/3), and 
obtain our final formula 



<f>p(x,y,z) = 8nu x u y —3 cos(2nu x px) cos(2nu y qy) + 

p,q>o Ppi 

4iru x u y Y —3 — cos(2nu y qy) + Anu x u y ^ — — cos(2-Ku x px) - ^ > 

q>0 Pi p>0 Pp 

2rru x u y \z\ + 2iru x u y P~ 1 + ® Q (P) ■ 



It has a singularity of 2nu x u y j3~ 1 which is independent of the particle coordi- 
nates. However, once the sum of (pp is taken over all particles, the singularity 
vanishes via the charge neutrality condition. For the other parts of Eq. (8) 
taking the limit (3 — > is trivial. The resulting formula (see Sec. 2.3) can be 
evaluated linearly in time, just like the Fourier part of the Ewald sum. This 
will be shown in more detail in section 4. 

Unfortunately the convergence becomes slower with decreasing z and for z = 
the sum is not defined. Thus we will need an alternative method for small z. 



2.2 Transformation of 4>p for z m 



For small particle distances, the term for k = I = is dominant (and must 
be omitted for the interaction of a particle with its own images). Therefore 
we leave it out for now and concentrate on a rapidly convergent formula for 
<j>p(x, y, z). This requires a little more work. For a more detailed derivation see 
[22,25]. 

Since we omit the k — I — term, the area to sum over has a hole. To 
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efficiently treat the sum over this area we split <f>p(x, y, z) — Ei + S 2 where 



Si = EE— and s * = £ 



2^0 fcez 



r k i 



Graphically this can be displayed like that: 



H 



We start by calculating Ei. Using the same argument as for formula (8) we 
obtain 



P -Pr k i 

Si = E E = 2«« E E MMe 2niu * px 

z^ofcez r ' kl 1^=0 P ez 

= 2u x E K (Me 2,rto + 2^E K o(^)- 
«,p^0 z^o 



While the first sum converges fast, the second one is still singular in (3 and 
has to be investigated further: 

2^E K o(« = 2u x J2Ml3pi) -2u x K {(3p) 

Z^O ZeZ 

= 2u x E K (/3p,) - 2u x (\og 2 - 7 - log(/3p)) + O {(3 2 ) 

where we have used the asymptotic behavior K (x) = — log | — 7 + C(a; 2 ) for 
x -> 0. 

The first sum can now be Fourier transformed again: 



E MP Pi) = ™ y E —5— (e 2mUyqy + e" 2 ™^) + tt % - 



-01*1 



<j>0 



,-Znu y q\. 



= 2nu v Re ( - 



2niu y qy 



nuytf- 1 - \z\) + O q {(5) 
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If we set C : = 2ixu y {\z\ + iy), we obtain 

-2nu y q\z\ -q£ ( L 

2™, E ^ z 2muyqv = E — = - MO + 1 - E ^tC" 



where 6 n are the Bernoulli numbers. This series expansion is valid only for 

2 



\(\ < 2n, which is fulfilled if \z\ < The last equality is found by integration 



from 



/ OO L \ I OO t 

1 E^ n -i =~ + E^ 2B - : 



Vn=0 



n! ; 2 ^(2n)! 



where the defining series -^—^ = S^Lo ^t^™ f° r the Bernoulli numbers was 
used. 

Using Re(— logC) = — log |C| = — log(27ru y p) and Re (J^J = nu y \z\ we obtain 
E MPPi) = - E -^T Re (C n ) - log(27r M ,p) + nuyp- 1 + (0) . 

l& n>2 nnl ^° 



It is easy to see that ( can be replaced by £ := 2nu y (z + iy) without changing 
the value of the sum. This is of advantage for the calculation of the forces by 
differentiation. 

Combining everything we obtain 

Si = 2u x E K ((3 pPl )e 2 ™* px -2u x J2-^i Re((2iru y (z + iy)) n ) - 

l,P^0 n>2 nn] (9) 

2u x \og{A-Ku y ) + 2u x log(/3) + 2-Ku x u y f3- 1 + 2u xl + ^(/J) . 



Now we concentrate on S 2 = Dfc>o ^ — l~ Dfc<o 6 rfc ■ It is sufficient to inves- 
tigate the first of the two sums, because by replacing x by —x the value of the 
second sum can be obtained. 

We start with 



E e fc e 10 f ^ ^ 1 E e n 

k>0 Tk k>0 ^ rfc k>0 



k>0 ^ rfc k>0 

/ 1 1 \ e-/ 3 '**; 
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Details about the precise derivation of the equation can again be found in 
[22,25]. 



Moreover by r k — k\ x = O (1) and log (1 — e a ) = logo; + 0(a) we obtain 

k— >oo 



To evaluate the sum J2k>o — ~kir) > we cons ider a TVy, such that > u x p + 
1. Then for k > we have p/|rr + fcA^] < 1 and therefore we can use the 



binomial series for < / 1 + , , 2 to obtain 



fc>AT^ 



fc^JVy, \ l X + kA x | ^x/ n>0 \ n J k>N^ \ X + k\ x \ 2n+1 

,(0),AT \ ^ /-|W (2n) (^+ M;c x) . , 2n 



n>0 



n J X 2 x n (2n)\ 



where ip^ are the polygamma functions. For details on these functions, see 
[24]. 

In summary we obtain for S 2 : 



(^ {2n) (^ + u x x) + ^ 2n \N^ - u x x)) 
V n J ( 2n ) ! 

2u x 7 + 2u x log(u x )-2 U:e log(/3) + X! - + +Pn^ lo S^))- 



Combining the formulas for S x and £ 2 we obtain, for |x| < 4f , |y|, |z| < 4f 
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and such that > u x p + 1, the final formula 

<f>/3(x,y,z) = Au x (Ko(pppi) + K (^pp_/)) cos(2nu x px) - 

l,p>0 

(^ 2 ")(A^ + Ux x) + ^ 2n \N^ - u x x)) 



n>0 



n J (2n)\ 



2u x \og(An^) + 2iru x u y [3- 1 + O (PlogP) . 

U x p^O 



(10) 



Of course this formula leads to the same singularity in (3 as formula (8) and the 
charge neutrality argument also holds for any combination of the two formulas 
as long as the sum is performed over all particles. 



2.3 Energy expressions 

Now we want to give the full expressions for the energy after taking the limit 
P -> 0. 

Setting 

fpq = \J(u x p) 2 + (u y q) 2 , f p = u x p , f q = u x q , 
Up = 2iru x p and u q = 2iru y q 

we obtain for z ^ 

e -2irf pq \z\ 

<f>(x, y, z) = Au x u y 7 cos(u p x) cos(uj q y) + 

p,q>0 Jpq 

( e -^h\A e -^U\A \ 
2u x u y cos(u q y) + Y 7 cos(^ p :r) - 2nu x u y \z\ . (11) 

\g>0 JQ p>0 JP ) 

Eq. (11) will be called the jar formula, since it will be used only for particles 
that are separated in the z-plane by a fixed minimum distance. See Sec. 4 for 
details of the implementation. 



10 



For \z\ < -f- we have 



(f>(x,y,z) =Au x (Ko(^pi) + K (cj p p_/)) cos(u p x) - 



l,p>0 




2n 



) 



2^l0gf47T^ . 



(12) 



In the same spirit as above this expression will be called the near formula. 

The far formula looks the same as an expression derived by Smith [15] and 
which is attributed to an approach originally due to Nijboer and de Wette [13]. 
Smith employed a spherical limit instead of a convergence factor, and therefore 
obtained a totally different singularity. Therefore his formula can be combined 
with any other formula using a spherical limit for small z, for example 2D- 
Ewald [11]. This combination leads then to the algorithm described in detail 
in Ref. [15]. Compared to this method, our combination has the advantage that 
the near formula is faster to evaluate than the 2D-Ewald sum and possesses 
an easy to use error formula that will be derived in the next section. Therefore 
the MMM2D algorithm is faster than the Smith algorithm, and in addition 
possesses full error control. 



2.4 The force 



Since the sums in equations (11) rsp. (12) converge absolutely, the electrostatic 
force Fi = —V Pi Ei can be derived by simple term- wise differentiation and the 
force can be calculated as 



N 
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where for z ^ F — (F x , F y ,F z ) T is given by 



e -2nf pq \z\ 

F x (x,y,z) = 8nu 2 x u y ^ V 7 sin(w p a;) cos{y> q y) + 

p,q>0 J PI 

Anu x u y ^2 e ~ 27r/p|2 ' sin(uj p x) , 

p>0 



e -2wf pq \z\ 

F y (x,y,z) = S-KU x u 2 y Y Q 7 cos(u p x) sm(u q y) + 

p,q>0 Jpi 

4nu x u y e~ 27r/ « |z| am(u q y) , 

q>0 



F z (x, y, z) = 8n sign(z)u x u y e 2n fp*\ z \ cos(uj p x) cos(uj q y) + 

p,<?>o 

An sign(z)u x u y ^ e~ 2nfq \ z \ cos(uj q y) + 

q>0 

Aix sign(z)u x u y ^2 e~ 2n f p ^ cos(uj p x) + 2ir sign(z)u x u y 

p>0 



For \z\ < 4f we have 



— ^— (x,y,^) T (s,y,z)^ (0,0,0) 
F(x,y,z) = F(x,y,z) + { (* 2 +2/ 2 +2 2 )s 

^0 (x,j/,z) = (0,0,0) 



where 



F x (x,y,z) = 8ttu x ^ p(K (w p pz) + K (o;pp_/)) sm(w p x) + 

Z,p>0 

iVy,-l / \ 



k=l \ r k 



1 \ ty(2n+l) {N + ^ _ ^(2n+l) (Ar _ UxX ^ 

S I » J w Kp) 



n>0 
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F y (x, y, z) = 8nul V p( + fa "^ )Kl( ^' i) 



cos(wpx) 



N, h -1 



E (iyy MOW* + ^.)) 2 "" 1 ) + E ( S + J- ) + 

1 (N^ - u x x)) 



3 x - (-\\ (V> (2w) (^ + ^x) + 4 2r - 
^^{n) (2n-l)! 



2(n-l) 



F z (x,y,z) =8iru 2 x ^ p /f^X^pP*) + ^ K i(^ P P O A cos ( UpX ) _ 
i,p>0 ^ ^ p ~ l ' 



4 ™^* E (g)T Re((2^(, + iy)) 2 *" 1 ) + E ( ^ + ^7 ) + 



n>l 



fc=l 



UxZ ^A n ) (2n-l)! 



(7V^ - u^)) 



(«xP) 



2(n-l) 



where K 1 is the modified Bessel function of order one. 



3 Error Estimates 



For an implementation which is of practical use, error estimates are needed. 
Since we calculate the energy (forces) by summing up pairwise energies Ei 
(their differential), it is reasonable to derive an upper bound for the error of 
E. The maximal pairwise error €e for the energy is induced by the calculation 
of the formulas (11) and (12) with finite cutoffs. Likewise we can derive an 
maximal pairwise error €f for the forces by formulas (13) and (14), and one 
can furthermore give an upper bound for the commonly used RMS-error of the 
forces. We will show later that the error distribution for MMM2D is highly 
non-uniform, and leads typically to an RMS-error is much lower than our 
bound. Thus the RMS-error is not the optimal error measure. Furthermore 
our implementation shows that an increase in precision has little impact on 
the calculation time of MMM2D, quite contrary to mesh based [8] and other 
methods[10]. 



3. 1 Error of the far formula 



For the far formula given by (11) and (13), respectively, we use a radial cutoff. 
So the summation is not performed over all (p, q) ^ 0, but only for those 
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(p, q) E T R where 



= {(p, q) e [0, oo) 2 \{0} | u 2 x (p - l) 2 + u 2 y (q - l) 2 < R 2 } . (15) 



By a simple approximation of the sums by integrals it is easy to derive upper 
bounds. For the potential we find the estimate 

(16) 

This upper bound is not valid for the forces. For all three components F x , F y 
and F z of the forces and the potential we find the upper bound 

e -2*R\z\ ( 1 \ 

t f = 2-kR + 2{u x + u y ) + — . (17) 



This value is precise only for F z . The other force components and the poten- 
tial show even a better convergence. Common to both error formulas is the 
exponential decrease of the error with R\z\. 



3.2 Error of the near formula 



The near formula given by equations (12) rsp. (14) contains three sums with 
different cutoffs. 

For the first sum containing Bessel functions it is reasonable to sum over all 
(j), q) G VLl where 



n L :=l(p,l) 



0<p<— und 0</< — + 1 1. (18) 



7TU X LU p 



If L > max(3uy,7iu x + u y ), \x\ < X x /2, \y\, \z\ < X y /2, the inequality Kq{x) < 
e~ x for x > 3 can be used to derive the upper bound 



/ gTiiAj / L 4- u \ ^ n " x 1 ^ 

r B = 8u x m&x(2nu x , l)e~ x « L y - 1 + V pe ^ u ^ p 

\ nu x X v \ ttu x J ±-1 



x /\y \ I u, x / 



(19) 



again by an approximation of the sum by an integral. The Gaussian bracket 
[•] denotes rounding up. This is an estimate for the absolute error of all force 
components and the potential. 
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For the second sum in (12) containing the Bernoulli numbers we use the esti- 
mate |&2n| < 4 (2^)2" fr° m [24] to obtain the overall error estimate 

r K = \Qu x Uy(u y p) 2N - 1 < 16V2u x u y 2~ N . (20) 

for the summation of this sum up to N. Note that this sum does not contribute 
to F x , due to the artificially broken symmetry in the derivation. The error 
estimate contains the particle position dependent term p. By using a table 
lookup scheme for N(p) one can chose the appropriate cutoff at runtime to 
speed up the calculation. 

For the last sum containing the polygammma functions there is no error es- 
timate necessary. The terms in this sum drop monotonously with alternating 
signs and therefore the absolute value of the last added term is an upper bound 
on the error. So the cutoff is determined at runtime. 



4 Efficient Implementation 

In this section we want to describe how the formulas (11) and (12) can be 
efficiently implemented to achieve the time scaling of 0(N^ log(iV) 2 ). 

The simulation box is split into B equally sized slices along the z-axis. We will 
see later how the parameter B should be chosen. Slice S contains all particles 
% that fulfill Jr < Zi ~* min < ^±1, where z m i n and z max are respectively the 
minimal and maximal ^-coordinates. Furthermore we set \ z = z max — z min 
and Ig should be an arbitrary enumeration of the particles in slice S. 

Now for all particles in slice S the interaction (potential or force) with the 
particles in the slices S — 1, S and S + 1 (if existent) will be calculated using 
the near formula (12) rsp. (14). For the other slices the far formula (11) rsp. 
(13) is used. In the following picture the splitting of the calculation for all 
particles in the black slice is shown: 




For the the near formula to be valid we need \z{ — Zj\ < \ y /2 for particles i 
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and j located in adjacent slices. This gives the constraint 



2A 2 A,, . . 

~~b £ f (21) 

Therefore \ y should be as large as possible and it might be useful to exchange 
the x- and y-axis if \ x » X y . In typical cases this will not be necessary since 
B is normally chosen large enough to yield a minimal computation time. 

The minimal distance of two particles that are treated by the far formula is 
b = X z /B. With this distance it is easy to use the error estimate (16) rsp. (17) 
for the far formula to calculate the cutoff R. Essentially we have 



4-1 Partial ordering of the particles 



First the sets I s are determined. This can be done in time O(N), as the 
following algorithm shows: 

Program 4.1 (Partial sorting) 

For all particles j determine its slice Sj as 

g q Z% Z-min 

Z-max Zmin- 

Add j to ■ 

The Gaussian bracket |_-J denotes rounding down. It is more efficient to reorder 
the particle list so that the particles of a given slice are stored adjacently, since 
later many indirect accesses via the sets Is can be omitted. 



4-2 Implementation of the far formula 



The calculation of the far formula consists of summing terms with frequencies 
(j), q) in the Fourier space. Now we show how the calculation of the contri- 
bution to the energies by a fixed frequency (p, q) can be done such that 
the computing time is O(N). The algorithm presented can easily be adapted 
for the calculation of the sums with only a single Fourier frequency and the 
|z|-sum. The calculation of the forces can be done simultaneously with the 
calculation of the potential using the same expressions. 
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In the beginning we concentrate on a single particle % located in slice Si. The 
far formula is used to calculate the contributions from all particles in non- 
adjacent slices, i. e. slices Sj where \Sj — Si\ > 1. First we restrict attention 
to the particles j in slices Sj < Si — 1, that is, to the particles in the set 
Us<Si-i Is- The cosine terms are separable via the addition theorem and since 
\zi — Zj\ = Zi — Zj for these slices, we have 

QiQj f cos{uj p (xi - Xj)) cos{uj q (yi - yj)) = 

jeis * pi 

S<S;-1 

qi — 7 cos(uj p Xi) cos(uj q yi) qje 27TfpqZj cos(u p Xj) cos(uj q yj) + 

I pi jei s 

S<S t -l 

^~^ 7V fpq z i 

qi — cos(ujpXi) sin(uj q yi) qje 2nfpqZj cos(u p Xj) sm(u q yj) + ^3) 

I Pi jeis 1 ^ 

S<S t ~l 

qi — 7 sm(u p Xi) cos(u q yi) ^ qje 2nfpqZ] sin(u p Xj) cos(uj q yj) + 

Jpq jeIjS 

S<Si-l 

qi — 7 sm(u p Xi) sin(u q yi) qje 27TfpqZj sm(uj p Xj) sin(uj q yj) . 

Jpq jeIs 

S<Si-l 



The sum for the S > Si + 1 is very similar, just the sign of the Zi and Zj is 
exchanged. For all particles j only the eight terms 

^(±,s/c,s/c) _ ^ e ±2n f pq Zj g j n j C08 ( UpX: .} gi n j C os(uJ q yj) 

are needed. The upper index describes the sign of the exponential term and 
whether sine or cosine is used for Xj and yj in the obvious way. These terms 
can be used for all expressions on the right hand side of Eq. (23). Moreover it 
is easy to see from the addition theorem for the sine function that these terms 
also can be used to calculate the force information up to simple prefactors that 
depend only on p and q. Once these terms are calculated for every particle, 
the rest of the calculation only consists of adding products of these terms 
^{±,s/c,s/c) w ^ correc t prefactors to the potential and the force. While it is 
easy to find these prefactors from the above formulas, it requires some care to 
insert the correct signs from the addition theorems into an implementation. 
Again we present the algorithm in pseudo code: 

Program 4.2 (Calculation of the far formula for a fixed (p, q)) 

(1) For all particles j calculate ^. ± ' s/ ' c ' s/ '^ '. 
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(2) For all slices I calculate 

~(±,s/c,s/c) _ \p (±,s/c,s/c) 
X\ / , Xj 

(3) For all slices k calculate 

A+,s/c,s/c) 
J- ,s / c,s / c) 

(4) For all particles i add the appropriate linear combination of the terms 
^Xt,s/c,s/c) ^ an( ^ f erms j an d foe terms ^' s ^ c ' s ^ c ' 1 (right hand terms) to 
Ei and Fj. Remember that a term with upper index + is always combined 
with a term of upper index — and vice versa. 

There are some tricks to speed up this algorithm. An easy way is to pre- 
calculate the sine and cosine terms for several frequencies p and q. And now it 
becomes clear that these loops are more efficient, as mentioned above, when 
particle data and results are stored consecutively in their respective arrays. 

If larger cutoffs R are needed, the terms e^ pqZi may get large rapidly. Then 
it might be useful to use the scaled coordinates e^ pq( - Zi ~ Zb ^ instead, where Zb 
is constant for a slice (e.g. the upper, lower border or the center). Then the 
terms x f° r the slices are calculated as before, the scaling is corrected during 
the summation for the £. The advantage of this is that the expression e^ pq ^ Zb ~ z ^ 
is numerically much more stable than e^ pqZb e~^ pqZ t> . 

Remember that this algorithm has to be executed for every frequency (p, q) 
and also for the single frequency sums and for the |z|-sum. 

4-3 Calculation of the near formula 

The calculation is algorithmically easy: It consists of two nested loops: 

Program 4.3 (Calculation of the near formula) 

For every slice S let Is = . . . ,ii} . Then for — ii, . . . , i\ add the self 
energy, i. e. the potential, that it feels from its own copies qf 0(0, 0, 0). 

Furthermore calculate the potential and the force the particle ik feels 
from the particles j = ik+i, ■ ■ ■ ,ii and the particles j G Is+i- Add these 
terms to the potential and the force of both particles, changing the sign 
of the force for particle j . 



jei s l<k-i 
S<k-1 

^(-,s/c,s/c) _ ^ ^(~,s/c,s/c) 

jeis l>k+l 
S>k+1 
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The self-energy can be calculated by formula (12). For X x = \ y = 1 it is 
0(0,0,0)^-3.90026. 

Now we want to give some hints for the calculation of the somewhat unusual 
functions involved. For the Bessel functions one can find approximations for 
example in [24]. Code using high order Chebychev approximations to calculate 
the Bessel functions is freely available, for example from the GNU scientific 
library [26]. 

Also we need a table of the Bernoulli numbers, which is easy to obtain. But 
it is more efficient to store the terms 

hn (2tt) 2 ". 



2n(2n 



The calculation of the sum over the Bernoulli numbers can be accelerated by 
making the cutoff depend on u y p using the error formula (20). 

The polygamma functions needed for the last sum are not calculated directly, 
instead we will give quickly converging power series for the functions 

i, (2n \x) := (^ 2n \N + x)+ ^ 2n \N - x)) and (24) 

^ 2n+l \x) := ^ ^ 2n+1 \N + x)- ij (2n+1 \N - x)) , (25) 

since these functions have a much better convergent series expansion than the 
polygamma functions. 

Letting 



y 2 W£\x) and ^\x)=( ^ /2 X +1) (x), 



the terms needed in the last sum can be calculated as 

u x ^ 2n \x)(u x p) 2n , u 2 Ji 2n+1 \x)(u x p) 2n , 
2nyu 3 Ji 2n) (x) (u xP ) 2{n ~ l) , 2nzu 3 x fy 2n) (x) (u x p) 2 ^ 



For the calculation of the power series we need the generalized zeta function 
or Hurwitz zeta function 

oo 

C(m,N) :=£r m . 

l=N 

Algorithms to calculate this function can be found in [24], but code calculating 
this function is also freely available, again from the GNU scientific library [26]. 
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Again from [24] we obtain with a little effort the series expansions 
ffi\x) = -2 ± ( 2n + ^ C(2n + 2k + 1, N)x* k , 

oo 

A 0) (x) = - 7 - 2 £ C(2A; + 1, iV)x 2fc , (26) 
fc=i 

4 2n+1) (a;) = 2* £ ( 2n + 2 + \ + X ) C(2n + 2k + 2, A^ . 
^.^ Computation time 

Using the implementation scheme given above, an estimation of the asymptotic 
computation time can be derived. We assume that the particles are distributed 
uniformly in a box. In other words, we assume that N/ B particles are located 
in every slice. 

The program 4.1 has clearly a computation time of O(N). 

The program 4.2 is executed for all frequencies (p, q) G Tr, that is 0(R 2 ) 
times. The first two steps of the program obviously have a linear computation 
time, the third step a computation time of order 0(B 2 ). The last step is linear 
again. In total we obtain a time needed for the program of 0(R 2 N)+0(R 2 B 2 ). 

The time for the calculation of the near formula for a single particle pair is 
assumed to be constant. This is sensible since at least the time for the sum 
over the Bessel functions is constant and the times for the other two sums are 
bounded. Therefore the time needed for program 4.3 is 0(N(2N / B)). 

If we insert R ~ B/\ z \og(B/(e A 2 )) from equation (22), we obtain an asymp- 
totical computation time of 

0[B 2 \og[B) 2 N) + 0(B 4 \og(B) 2 ) + 0{N 2 /B) + 0{N) . (27) 

for the complete program. 

Choosing B ~ A^, we obtain a computation time of 

G(nI log(A0 2 ) + 0{m \og(N) 2 ) + 0(Nl) . 

This gives a total computational time of 0(N§ \og(N) 2 ) : which is easily seen 
to be optimal from equation (27). 

Although we know the optimal proportionality of B, this parameter has to be 
tuned to the underlying hardware. But this can be done by a simple trial and 
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error algorithm at runtime, since the computational error does not depend on 
B. 



5 Numerical Examples 

In this section we give some numerical results from our implementation of 
the MMM2D algorithm. The computations were performed on a Compaq 
XP1000 with a 667 MHz EV67-CPU. 

5. 1 Regular Systems 

First we want to demonstrate that MMM2D reproduces some known results 
from Ref. [10]. We also give the consumed computation times for MMM2D 
and for an implementation of the 2<i-Ewald method. 

We consider the following systems, which are in more detail described in 
Ref. [10]: 

System 1: 100 particles distributed regularly in a square of length L — 1. The 
particles have a unit charge value with chessboard-like alternating signs. 

System 2: The same as above, but one particle elevated by 0.5 out of the plane. 

System 3: This system consists of 25 particles in a square and one elevated by 
0.2 out of the plane in the center of the square. The charges in the plane are 
again chessboard-like distributed, the full system is charge neutral. 

Table 1 shows some results for these systems. For the numerical values, we 
tuned MMM2D for a maximal pairwise error of e = 10~ 6 in the forces by 
the error formulas given above. For the timings Tmmm2D we used a maximal 
pairwise error of e = 10~ 4 (timings for e = 10~ 6 in parentheses). The 2<i-Ewald 
method was tuned for a RMS-error of 10~ 4 , since for a typical simulation an 
error of 10~ 4 is sufficient. Note that for MMM2D the actual RMS-errors 
were 10~ 7 for e = 10 -4 and 10~ 9 for e = 10 -6 . For system 1 the forces on 
all particles are equal in magnitude and the force on an arbitrary particle is 
given, for Systems 2 and 3 the forces on the elevated particle are given. 

In system 1 the z-distance is for every pair of particles, and therefore only 
the near formula is used. For the systems 2 and 3 only the far formula was 
used for the interactions with the elevated particle, and the precision for those 
systems is much better than for the first, as can be inspected in Tab. 1. This is 
due to the exponential factor containing the z-distance in the error formulas 
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Table 1 

Force and energy results and computation times for the systems 1, 2 and 3 



System 


F x 


Fy 


F z 


Etotal 


TMMM2D 

[ms] 


T2d-Ewald 

[ms] 


1 


5 


1(T 9 


2 • i(r 8 





-807.771313 


74.1(110) 


318 


2 


-2 


■ io- 17 


-2 ■ 10~ 17 


-7.765381 


-792.588065 


72(108) 


317 


3 


4- 


1Q -16 


1 • 1(T 15 


-10.364160 


-86.565859 


4.91(7.32) 


21.1 



Eqs. (17) and (16). This factor leads to a non-uniform error distribution with 
respect to z and forces MMM2D to be extremely precise for some particle 
positions, and therefore the forces of the elevated particles are calculated with 
likewise precision. We will show this in more detail in the next part. The 
computation times will be discussed later. 

5.2 2-particle systems 

For Figure 1 we used two particles, one located at (0,0,0.5), the other ran- 
domly placed in a cubic box of unit length. We show the error in the z-force 
on the fixed particle for a calculation with B = 8 and a maximal pairwise 
error of e = 10~ 6 . 

o.oooi 

le-05 

le-06 

fc N le-07 
< 

S le-08 
t3 

u le-09 

u 

<2 le-10 
le-11 
le-12 
le-13 

0.2 0.4 0.6 0.8 1 

z 

Fig. 1. Absolute z-force error in dependence of the z-distance of the two-particle 
system (compare text) 

Since the box has a height of 1 and B = 8, the slices have a height of 0.125 
and the forces are calculated with the near formula if and only if 0.375 <= 
z < 0.750. In this distance range is the maximum error nearly constant and 
scatters well below the maximal pairwise error. The additional structure which 
is observable in the error distribution is due to the way the different sums in 
the near formula (12) are truncated depending on z in order to reach the 
prescribed precision. 
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For z < 0.375 the exponential drop of the error predicted by the far error 
formulas Eqs.(17) and (16) is clearly visible. This effect cannot be seen on the 
right side for large z since the error of the far formula is symmetrical, while 
the distance where the near formula is used is asymmetrical with respect to 
the fixed particle. If the fixed particle would be located at the upper bound 
of its slice, for example at (0,0,0.624), the far error would be visible on the 
right side, but not on the left. 

The exponential drop of the far formula error suggests that one could reduce 
the numbers of added Fourier components for more distant slices to reduce 
the computation time. However, this would only save a number of operations 
(additions and multiplications) of the order of the number of slices, and the 
overhead for singling out the slices compensates the gain in time. 

In the time scale considerations above, the error estimate for the far formula 
has a central place. Therefore we want to show that the estimates Tp (for 
the forces) and te (for the energy) are accurate. Figure 2 shows the error 
estimates and the maximal errors in dependence of the cutoff R. The system 
investigated here again consists of two particles, one fixed at (0,0,0), the 
other one randomly placed in [0, 1] x [0, 1] x [0.25, 1]. It can be seen that both 
error formulas overestimate slightly the error, but are sufficiently close to the 
numerical results. 

l 

0.01 
0.0001 

a 

§ le-06 
I 

o le-08 

■a 

le-10 
le-12 
le-14 

2 4 6 8 10 12 14 16 18 20 
R 

Fig. 2. Maximal absolute errors of the z-force and the energy for the 2-particle 
system (compare text) 

5.3 Randomly distributed systems 

Now we want to investigate the RMS-error scaling of MMM2D and the 
computation time for larger systems. The systems investigated consist of iV 
particles placed randomly in a cubic simulation box of unit length. Half of the 
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particles carry a charge of 1, half of them a charge of —1. For the calculations 
MMM2D was tuned to a maximal pairwise error of e = 10~ 4 . 

The error distribution is highly non-uniform as we have shown before. There- 
fore it is more safe to use the error formulas given above to constrain the maxi- 
mal pairwise error rather than the RMS-error, as is done usually. Nevertheless 
we investigated the actual RMS-error for the uniformly random systems de- 
scribed above. Figure 3 shows the RMS-error for different particle numbers 
and slice numbers. By the arguments given in Sec. 5.2, the error drops when 
more slices are used. 



B=4 




10 100 1000 

N 

Fig. 3. RMS-error of the forces calculated by MMM2D for different values of B 
in dependence of the particle numbers N for fixed e = 10~ 4 

If one assumes a standard deviation of the pairwise error of a, and furthermore 
that the pairwise errors are independent, which is true for random systems, 
it is easy to see that the RMS-error has an expectation value of [27], 
with Q 2 = Y.i if, which can be seen in Figure 3. Thus our RMS-error behaves 
like the ones for the usual Ewald methods [27]. Furthermore we have a < e 
so that the RMS-error can be bounded using our error estimates. Note that 
again the RMS-error is much smaller than the maximal pairwise error (for 100 
particles, 3 • 10~ 6 instead of 10~ 4 ), although the gain is smaller than for the 
highly structured systems 1-3. 

Having laid the theoretical foundations of the time scaling in Sec. 4, we want 
to show that the theoretical scaling can be achieved in a real computation. 
Figure 4 shows the optimal computation time depending on the particle num- 
ber N for the systems considered above. We also give data for different box 
heights Az. The x and y sides are still fixed to unit length. The pairwise er- 
ror in the forces was fixed to e = 10~ 4 , resulting in the better RMS-errors 
given above. The optimal number of slices B depends on Az. For example, 
for Az = we always have B = 1 and a purely 0(N 2 ) algorithm, as is seen 
in the figure. For Az = 1 the optimal slice numbers were B = 4 for N < 100, 
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B = 8 for N < 1000 and B = 16 for N = 4000, and the 0(N 5 / 3 ) scaling 
is achieved. It can also be seen that MMM2D becomes faster with growing 
Az, as one expects from the error formula. For comparison also computation 
times for the 2o?-Ewald method (at Az = 1) and the theoretical scaling curves 
of 0(N 5 / 3 ) and 0{N 2 ) are included. 

The computation times show that MMM2D is always faster than the HBC 
method, even for small particle numbers (although we achieve a much smaller 
RMS-error). This can also be seen in the computation times given in Tab.l. 
These computation times also demonstrate that increasing the precision does 
not massively impact the computation time of MMM2D. This is due to the 
exponential error decrease with respect to the cutoff R. Since System 1 in this 
table is purely two dimensional, MMM2D only uses the near formula for 
the calculation. Therefore even our near formula is superior to the 2<i-Ewald 
method. 

Since MMM2D and the algorithm described by Smith [15] are similar up 
to the difference that Smith uses the 2o?-Ewald method instead of our near 
formula, MMM2D is also faster than the algorithm of Smith. From the num- 
bers given by Widmann [10] one can furthermore conclude that MMM2D is 
also faster than the algorithm proposed by Hautman and Klein[12], although 
we did not test this method explicitly. 

MMM2D is also similar in speed at our accuracy with respect to the two 
recent Ewald-type algorithms proposed in Refs. [16,17]. However, these two 
methods also have the drawback that no error estimates exist, and these will 
also be hard to develop because of the numerical integration involved. Fur- 
thermore it seems unreasonable that their algorithm gets slower for larger Az, 
while MMM2D gets faster, and that their algorithm is fastest for the largest 
tested real-space cutoffs[16] although they claim that their new Fourier-space 
expression is responsible for the increase in speed. On the basis of their num- 
bers and non-comparable test cases we are unable to draw further conclusions. 



6 Conclusion 

We developed a new method termed MMM2D to accurately calculate the 
electrostatic energy and forces on charges being distributed in a two dimen- 
sional periodic array of finite thickness (2d+h geometry). It has the advantage 
of requiring only a computational time of C(iV 5//3 ) instead of a quadratic scal- 
ing as previous methods show. It is based on a convergence factor approach 
and as such does not require any fine-tuning of any Ewald parameters for con- 
vergence. Moreover, we derived rigorous error bounds and numerically verified 
them on specific examples. 
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Fig. 4. Optimal computation times for different particle numbers N. 



Our method is not only asymptotically faster than all 2(i-Ewald methods, but 
is about four times faster than the HBC method for small numbers of parti- 
cles, and should also be superior to the method of Hautman and Klein, and to 
the new methods of Refs. [16,17] at our accuracies. Our energy and force re- 
sults achieve a maximal pairwise error of 10~ 4 with minimal cutoffs, and 10~ 8 
can easily be achieved without large increase in computational time. There- 
fore our method is even more superior whenever a high accuracy is needed. 
Although the presented formulas may look awkward at the beginning, once 
implemented, the method is easy to use because the cutoffs for the sums can 
be easily determined at runtime. Only the number of slices B must be op- 
timized for speed, but, although not implemented, this parameter could also 
be determined during execution. This is in contrast to basically all 2ci-Ewald 
methods, for which the tuning of the various parameters is very crucial and 
error estimates do not exist. 

The numerical code for MMM2D in C++ can be requested from the authors. 
This code also provides interfaces for use with FORTRAN and C. 
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A Equivalence of the sum with the convergence factor 



Let pel 3 and for k,l G Z let n k i '■— (k\ x ,l\ y ,0) T and n = \n k i\ = 
k 2 \ 2 x + l 2 \ 2 . Then it can be shown that for charge neutral systems 



o (V 3 ) 



\P + n kl \ \p-n kl \ n 
(see for example [5] or [25]). 

Now we rewrite the sum with the convergence factor as 

-0\pj+n k i\ 



N 



y^' Qj^_ 

5=0 k 2 +l 2 =S 3=1 \Pj + Ukl I 



N 



(A.1) 



(A.2) 



For any a n where /3|a„ — n| < 1 



„-/3(a„-n) 1 oo /_-iy 

^ — = ^-ESt-^K-")' 



(A.3) 



where 1 9((3, a n , n)\ < \ ^~^ , 



Therefore for sufficient small /3 we have 

e -/3|p+n fc; | g-/3|p-n fcJ | 2e~^ n 



1 1 



+ 



b + ra«| \p-n M \ 
n 



1 + 



+ e _/Jn \6((3, \p + n k i\,n) + 6((3, \p - n M \ , n) \ 



f3e 



-pn 



2 e -Pn 



0(1) + '—^0(1) + 



(3 2 e 



n 



n 



-Oil) 



This can be seen by applying equation (A.l) to the two first terms and equation 
(A.3) to the last. By approximating the sums over these terms separately as 
integrals it is easy to see that the sum in equation (A.2) converges uniformly 
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for small (3 and therefore 




\Pj + n ki \ \pj-n M \ n 



e -P\Pj+nu\ e -P\P]-n k i\ 2e _/3 '' 



\Pj + n M \ \pj-n M \ n 



1 1 2 



In 



) 



(A.4) 



oo JV 

= E E E' 



S=0fe2+j2=S j=l 



References 



[I] P. Ewald, Die Berechnung optischer und elektrostatischer Gitterpotentiale, 
Ann. Phys. 64 (1921) 253-287. 

[2] J. Kolafa, J. W. Perram, Cutoff errors in the ewald summation formulae for 
point charge systems, Molecular Simulation 9 (5) (1992) 351-68. 

[3] R. Sperb, An alternative to ewald sums - part i: Identities for sums, Molecular 
Simulation 20 (3) (1998) 179-200. 

[4] S. W. de Leeuw, J. W. Perram, E. R. Smith, Simulation of electrostatic systems 
in periodic boundary conditions, i. lattice sums and dielectric constants, Proc. 
R. Soc. Lond. A 373 (1980) 27-56. 

[5] S. W. de Leeuw, J. W. Perram, E. R. Smith, Simulation of electrostatic systems 
in periodic boundary conditions, ii. equivalence of boundary conditions, Proc. 
R. Soc. Lond. A 373 (1980) 57-66. 

[6] R. Sperb, R. Strebel, An alternative to ewald sums, part 3: Implementation and 
results, Molecular Simulation 27 (1) (2001) 61-74. 

[7] R. W. Hockney, J. W. Eastwood, Computer Simulation Using Particles, IOP, 
London, 1988. 

[8] M. Deserno, C. Holm, How to mesh up Ewald sums. i. a theoretical and 
numerical comparison of various particle mesh routines, J. Chem. Phys. 109 
(1998) 7678. 

[9] J. Lekner, Summation of coulomb fields in computer simulated disordered 
systems, Physica A 176 (1991) 485-498. 

[10] A. H. Widmann, D. B. Adolf, A comparison of ewald summation techniques for 
planar surfaces, Comp. Phys. Comm. 107 (1997) 167-186. 

[II] D. M. Heyes, M. Barber, J. H. R. Clarke, Molecular dynamics computer 
simulation of surface properties of crystalline potassium chloride, J. Chem. Soc. 
Faraday Trans. II 73 (1977) 1485. 



28 



[12] J. Hautman, M. L. Klein, An ewald summation method for planar surfaces and 
interfaces, Mol. Phys. 75 (1992) 379. 

[13] B. R. A. Nijboer, F. W. de Wette, On the calculation of lattice sums, Physica 
23 (1957) 309-321. 

[14] B. Nijboer, On the electrostatic potential due to a slab shaped and a semi- 
infinite nacl-type lattice, Physica A 125 (1984) 275-279. 

[15] E. R. Smith, Electrostatic potentials for thin layers, Mol. Phys. 65 (1988) 1089- 
1104. 

[16] M. Kawata, M. Mikami, Rapid calculation of two-dimensional ewald 
summation, Chem. Phys. Lett. 340 (2001) 157-164. 

[17] M. Kawata, U. Nagashima, Particle mesh ewald method for three-dimensional 
systems with two-dimensional periodicity, Chem. Phys. Lett. 340 (2001) 165- 
172. 

[18] A. Arnold, C. Holm, Chem. Phys. Lett, to appear. 
[19] J. Lekner, Physica A 157 (1989) 826. 

[20] R. Sperb, Extension and simple proof of lekner's summation formula for 
coulomb forces, Molecular Simulation 13 (1994) 189-193. 

[21] F. S. Csajka, C. Seidel, Strongly charged polyelectrolyte brushes: A molecular 
dynamics study, Macromolecules 33 (2000) 2728-2739. 

[22] R. Strebel, Pieces of software for the Coulombic m body problem, Dissertation 
13504, ETH Zuerich (1999). 

URL http : / / e-collection . ethbib . ethz . ch/ show?type=diss&nr=13504 

[23] F. Oberhettinger, Fouriertransforms of distributions and their inverses, Vol. 16 
of Probability and Mathematical Statistics, Academic Press, New York, 1973. 

[24] M. Abramowitz, I. Stegun, Handbook of mathematical functions, Dover 
Publications Inc., New York, 1970. 

[25] A. Arnold, Berechnung der elektrostatischen wechselwirkung in 2d + h 
periodischen systemen, Diploma thesis, Johannes Gutenberg-Universitat (may 
2001). URL 

http : / / www . mpip-mainz . mpg . de/ www/ theory/ phd_work/ dipl-arnold . ps . gz 

[26] GNU, Gnu scientific library (GSL) . 

URL ftp : //sources . redhat . com/pub/gsl 

[27] M. Deserno, C. Holm, How to mesh up Ewald sums. ii. an accurate error 
estimate for the p3m algorithm, J. Chem. Phys. 109 (1998) 7694. 



29 



